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Abstract 

In recent years, the availability of highly pure stable isotopes has made possible the investi- 
gation of the dependence of the physical properties of crystals, in particular semiconductors, on 
their isotopic composition. Following the investigation of the specific heat (Cp, C^) of monatomic 
crystals such as diamond, silicon, and germanium, similar investigations have been undertaken for 
the tetrahedral diatomic systems ZnO and GaN (wurtzite structure), for which the effect of the 
mass of the cation differs from that of the anion. In this article we present measurements for a 
semiconductor with rock salt structure, namely lead sulfide. Because of the large difference in the 
atomic mass of both constituents (Mpb= 207.21 and (Mg=32.06 a.m.u., for the natural isotopic 
abundance) the effects of varying the cation and that of the anion mass are very different for this 
canonical semiconductor. We compare the measured temperature dependence of Cp ~ C^, and the 
corresponding derivatives with respect to (Mpb and Mg), with ah initio calculations based on the 
lattice dynamics obtained from the local density approximation (LDA) electronic band structure. 
Quantitative deviations between theory and experiment are attributed to the absence of spin-orbit 
interaction in the ABINIT program used for the electronic band structure calculations. 

PACS numbers: 63.20.Dj, 65.40.Ba 
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I. INTRODUCTION 



Lead sulfide belongs to the family of the lead chalcogenides (PbS, PbSe, PbTe) which 
are small bandgap semiconductors that crystallize in the rock salt structure. They differ 
from the nowadays more commonly encountered tetrahedral semiconductors (e.g. GaAs, 
ZnO, with zincblende and wurtzite structure, respectively) in that they possess 10 valence 
electrons per formula unit, as opposed to 8 for the zincblende- and wurtzite-type materials. 
The lead chalcogenides are found as minerals in nature, bearing the names: galena (PbS), 
clausthalite (PbSe) and altaite (PbTe). Whereas galena is rather abundant and is found 
as large, highly perfect crystals, clausthalite and altaite are rare. Some semiconducting 
properties of PbS have been known since the mid 1800, in particular the rectifying prop- 
erties of metal-PbS contacts which were used as detectors in early radio receivers^. The 

n 

thermoelectric properties of PbS were reported as early as 1865 (Ref. |2). The photoresistive 
and photovoltaic properties were reported in 1930 (Ref. [Sj). A number of reviews which 
summarize the physical properties of PbS have appeared^i^. Particularly interesting is the 
anomalous dependence of the electronic energy gap on temperature^. 

The specific heat of PbS was first measured by Eastman and Rodebush above 63 K 
and subsequently by Parkinson and Quarrington for temperatures between 20 and 270K-i^. 
Such measurements were extended down to IK by Lykov and Chernik^. Although the typical 
carrier concentrations of the PbS samples were between 10^^ and 10^^ cm~'^, no evidence 
of a free carrier contribution to the specific heat was found, except possibly below 2K— . 
Thus the measured specific heat represents the contribution of the lattice vibrations vs. 
temperature, obtained at constant pressure (Cp). Most lattice dynamical calculations, such 
as those of Elcombe, represent the specific heat measured at constant volume (Cy)^. The 
difference between Cp and is given by^: 

Cp{T) - Cy{T) = al{T) ■ B ■ Vmoi ■ T, (1) 

where is the temperature dependent thermal expansion coefficient, B the bulk modulus 
and Vmoi the molar volume. According to Eq. (1), the difference Cp - C^ increases with 
increasing T. Replacing standard values of a;t,(T), B and Vmoi found in the literature we 
obtain from Eq. (1) the difference Cp - C^ = 0.18 J/mol K at the highest temperature used 
in our measurements (~280 K). This difference is smaller than typical error bars in our 
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measurements and will be neglected here. We shall therefore denote the specific heat, either 
theoretical or experimental, by ■ 



II. THEORETICAL DETAILS 

The specific heat of PbS vs. T has been calculated in the harmonic approximation from 
the vibrational free energy F{T) as obtained with the ABINIT code^^. This code derives 
the phonon dispersion relations from total energy calculations and uses them to evaluate 
F{T) with the expression: 



In Eq. ([2]), ks is the Boltzmann constant, ns the Bose-Einstein factor, and p{uj) the 
phonon density of states (PDOS). The high frequency cutoff of the latter defines the upper 
limit of integration in Eq. ([2]). 

The theoretical calculations were performed in the framework of the density functional 
theory (DFT) by using a local density approximation for the exchange-correlation Hamilto- 
nian as implemented in the ABINIT package. Kohn-Sham orbitals are expanded in plane 
waves with the use of the Hartwigsen-Goedecker-Hutter pseudopotentials to describe the 
valence electrons (4 electrons for Pb and 6 electrons for S)^^. A 70 Ry cutoff was used and 
a Monkhorst Pack grid of 8x8x8 points was used to describe the electronic properties in 
the Brillouin zone. Before vibrational calculations were performed, the cell parameter was 
optimized to a value of 5.808 A which was used thereafter. The ABINIT program determines 
the elements of the dynamical matrix by perturbation theory. They were calculated for a 
12xl2x 12 mesh and then Fourier interpolated to obtain the thermal properties discussed 
here. 

The ABINIT program determines the elements of the dynamical matrix by using pertur- 
bation theory. The various isotopic masses used in the calculation were introduced at the 
level of the dynamical matrix, assuming that the matrix elements of the latter do not depend 
on those masses. Possible dependence of dynamical matrix elements (i.e. force constants) 





From Eq. ([2]) the specific heat is obtained with: 




(3) 
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on mass can be estimated from Raman and inelastic neutron scattering measurements^. 
They are several orders of magnitude smaller than the corresponding isotope mass changes 
and thus can be neglected for the purpose of the present work. 

Under these assumptions and for monatomic crystals (diamond, Ge, Si) all frequencies 
Lo are proportional to M~°'^. Because of the appearance of the argument (uj/T) in ub (Eq. 
2), it is possible to relate the derivative of with respect to the isotopic mass M to that 
with respect to T. The following expression is found^ 

d\niC,,jT^) ^ 1 . dHC,,jT^) 

d\nM 2 ^ d\nT ^ ' 

In binary (or multinary) crystals, the masses of the different elements have a different ef- 
fect on each phonon frequency which, in general, cannot be evaluated with simple arguments. 
This effect depends on the phonon eigenvectors or, equivalently, on the DOS projected on the 
atom under consideration^. As we shall see below, for PbS because of the large difference 
in the masses of Pb and S, the lattice vibrations are predominantly Pb-like up to 120 cm~^ 
(acoustic phonons). Higher frequency vibrations (optic phonons) are overwhelmingly S-like. 
The frequency of the acoustic phonons thus should vary, to a good approximation, like the 
inverse square root of the Pb mass (Mpj^ ) whereas that of the optic phonons should vary 
like Mg-^/^ 

For computational reasons (the limited number of fc-points being sampled) the calculated 
specific heats are not reliable for T below ~5K. A similar low T limit also applies to the 
experimental data due to the small values of and the limited accuracy of the experimental 
procedure using sample masses of ~25mg for the isotope enriched samples. T < 5K corre- 
sponds to the Debye limit in which is proportional to . This expression enables 
us to derive for the mass derivatives in the limit T ^0 the following analytic expressions^^: 

dlnCjT^ _3 Mpb dlnCjT^ _3 Mg _ ^ 

dlnMpb 2Mpb + Ms ' rflnMs 2 Mpb + Ms ' 

The values given in Eqs. ([5]) can be used to extrapolate those obtained at 5K to T ^ 0. 
In order to achieve convergence of the heat capacity and its logarithmic derivatives with 
respect to Mpb and Ms a grid of 80x80x80 q points in the irreducible Brillouin zone was used 
for the evaluation of the PDOS and the thermodynamic properties. The number of q points 
used in the integrations (e.g., Eq. [2]) is particularly important for the convergence of the 
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calculations of the mass derivatives at low temperature. These calculations were performed 
by evaluating the differences of Cy for two different lead and sulfur pairs of isotopes. Below 
5K, the calculated derivatives depended strongly on the pair of isotopes and the sample 
chosen. We kept only the values for T as given by Eq. \5\ 

III. LATTICE DYNAMICS: DISPERSION RELATIONS AND DENSITY OF 
PHONON STATES 

Figure [1] displays the phonon dispersion relations of natural PbS obtained from the ab 
initio electronic band structure using the procedure described above (solid lines). Exper- 
imental points obtained by inelastic neutron scattering (INS)ii are denoted by diamonds. 
Of particular interest is the large LO - TO splitting which, by virtue of the Lyddane-Sachs- 
Teller relation^?, indicates a large value of the low-frequency dielectric constant (e(0) ~ 180, 
e(oo) ~ 18) which is related to the nearly ferroelectric character of PbTe^. 

The INS points in Fig. [1] agree reasonably well with the calculated dispersion relations 
except for the TO modes along the S and the A directions. Previous calculations were of 
a semi-empirical nature with parameters adjusted to fit the INS result s^ii^S. It is therefore 
not surprising that the agreement between the calculated dispersion relations and the INS 
results is somewhat better than that of Fig. [TJ 

We conjecture that the discrepancies between theory and experiments displayed in Fig. [1] 
may be due to the lack of the spin-orbit interaction contribution in the corresponding LDA 
Hamiltonian. This conjecture is supported by recent ab initio calculations of the phonon 
dispersion relations of bismuth performed with and without spin-orbit interaction^^, ab 
initio electronic band structure calculations which include spin-orbit coupling have been 
recently performed for PbS, PbSe, and PbTe but, to our knowledge, they have not been 
used to derive the phonon dispersion relations. Ab initio calculations of the LO and TO 
frequencies at q=0 have been reported for PbTe^^. 

The TO bands of Fig. [1] are nearly dispersionless along S and A. This fact suggests 
phonon modes in which the lead and the sulfur vibrations are only weakly coupled, a con- 
jecture which is confirmed by the calculated phonon eigenvectors and also by the PDOS 
displayed in Fig. [2] as projected on the lead ( predominant below 120 cm~^) and the sulfur 
(predominant above 120 cm~^) eigenvectors. The van Hove singularities of the TA modes at 
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the L and X points are responsible for the sharp, Pb-hke spike seen at ~50 cm~^ (cf. Fig. 
[2]), whereas the flat regions of the TO bands are responsible for the much lower peak found 
at ~130 cm~^. 



IV. EXPERIMENTAL PROCEDURE 
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. Two 



The preparation and some properties of the PbS samples are described in Ref. 
types of samples with the natural isotopic composition (Mpb =207.21 a.m.u., Ms =32.05 
a.m.u. ) were used. One was a piece of natural galena (from Creede, CO, USA, 16 mg) 
whereas the others were pieces of the same size taken from the synthetic ingots used in 
|. We also measured a piece of the '^'^^Pb^^S of Ref. 



Ref. 
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and a piece from an ingot 



prepared in the same way with the isotopic composition ^^^Pb'^'^^S. All samples were p-type 
(as determined from the sign of the thermoelectric power), the natural sample had a carrier 
concentration of 2x10^^ cm~^ whereas the synthetic ones had carrier concentrations around 

n 

10^^ cm ^ (as determined from the infrared plasma edge, see Ref. l26l ). 

The synthetic ingots were prepared by first reacting the corresponding pure elements and 
then subliming the compound in an argon atmosphere as described in detail in Ref. y. 

The heat capacities were measured between 5 and 280K with a PPMS system (Quantum 



Design, 6325 Lusk Boulevard, San Diego, CA.) as described in Ref. |l9|. All natural samples 
measured gave the same results within error bars. 

V. EXPERIMENTAL VS. THEORETICAL RESULTS 

Figure [3] shows the specific heat measured for our ^i^tpi^natg g^j^iples for T between 1.8 
and 270K and that reported in Ref. js for T between 20 and 260 K. The solid line represents 
the results of our ab initio calculation. Whereas the two sets of experimental data agree 
within error bars (the widths of the symbols in the plot), the calculated curve falls slightly 
below the measurements, especially around 80K. The sign of this deviation corresponds to 
calculated phonon frequencies lying above the experimental ones, as is the case for the TO 
modes in Fig. [H Notice that at temperatures above 280K Cy tends to the Petit and Dulong 
value {Cy ~ 49 Joule/mole K )^. 

In order to be able to read the low temperature data of Fig. [3l which according to Debye's 
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law should be proportional to T^, we display in Fig. 4 the quantity Cv(T)/T^ . This figure 
reveals that the region of validity of Debye's law is rather small (T < 3K). This follows from 
the fiattening of the TA dispersion relations with increasing q. Figure H] also shows that 
the calculated value of C^/T^ at its maximum {T^ax ~ UK) is about 23% lower than the 
measured one. This discrepancy can be assigned to that of the calculated TO frequencies 
(Fig. [1]) which has already been attributed to the lack of spin-orbit interaction in the band 
structure calculations. We recall that in the case of bismuth a similar discrepancy is found 
when the ah initio calculations are performed without spin-orbit coupling^. 

Using Debye's theorent^, Fig. Hlcan be recast into a form which exhibits the temperature 
dependence of an equivalent Od(T) (Fig. \5^. This figure, which includes calculated^i as well 
as experimental data, reveals a value of 6^ equal to 210 K in the T — > limit. Od{T) 
decreases to a minimum of 135 K at ~ 20K. 

Od{T) saturates at 9-q ~ 230 K for T > 230K. Note that in spite of the large scatter 
in our experimental points, the agreement between theory and experiment is rather good, 
except in the region around 75 K where a difference of about 1% appears. In order to get a 
feeling for the dependence of the results of Fig. H] on isotopic masses, especially around the 
maximum, we plot in Fig. [6^ and [Ub ah initio calculations for different sulfur (Fig. [6^) and 
lead (Fig. [6]d) isotopes. We note that a 2% variation in Mpb results in a 3% variation in 
C^/T^ whereas a 3% variation in Ms results in a 0.5% variation of C^/T^ . This reveals a low 
sulfur content in the eigenvalues of the acoustic modes as already discussed in connection 
with Fig. El 

Figure [7| displays the change in C^/T^ measured and calculated for a change of one a.m.u. 
in the mass of Pb in PbS. Below 6K, the calculations (and also the measurements) yield 
unreliable results. Correspondingly, the calculated curve has been extrapolated to the value 
calculated for T ^ using Eqs. [5] and the Debye limit for ( C^/T^ ~ 3/iJ/mol K^). The 
maximum in both, experimental and calculated data occurs for T ~ 12 K. Note that the 
larger isotope mass yields a larger value of C^/T^- 

Figure [HI shows experimental and calculated results similar to those in Fig. [7|but obtained 
for a change of 2 a.m.u. in Ms. The calculated curve displays a maximum at 35 K. The 
measured points scatter markedly around this temperature but, on the average, they are also 
compatible with the existence of a maximum at 35 K. Below ~10K both sets of data become 
unreliable. We have also extrapolated the calculated curve to join the value predicted with 
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Eqs. [5] and the Debye law. 



VI. DISCUSSION 

As shown in Figs. [3] and H] our calculations, based on phonon dispersion relations ob- 
tained from ab initio electronic band structures, reproduce reasonably well the experimental 
results. The discrepancy which appears at the maximum of Fig. H] can be attributed to the 
lack of spin-orbit interaction in the otherwise relativistic Hamiltonian used for the ab initio 
electronic calculations. The peak in the AC^/T^ of Fig. [7] occurs at 12 K, a temperature 
about one-third that of the peak seen in Fig. [HI The ratio of both temperatures approxi- 
mately equals the square root of the corresponding masses. This fact allows us to assign the 
peak in Fig. [7] to vibrations of lead atoms and that in Fig. [8] to those of sulfur atoms, in 
agreement with the assignment of the structures in the PDOS of Fig. [2] below and above 120 
cm~^. The Pb-like PDOS of Fig. [2] shows a sharp peak at 55 cm"^ (~ 80K in temperature 
units) which has been discussed in Sect. III. The S-like PDOS spreads over a broad band 
with a peak at 133 cm"^ (~ 196 K). The ratio of the temperatures of these two peaks to 
those of the maxima in Figs. [7^nd[8]is ~6. By modeling the two peaks in the PDOS with 
single harmonic oscillators (Einstein model) it is possible to calculate peaks in AC^/T^ at 
temperatures close to those displayed in Figs. [7] and [HI For simplicity, let us consider an 
Einstein oscillator with a frequency equivalent to lOOK which is inversely proportional to 
the square root of its mass. The logarithmic derivative of the corresponding specific heat 
with respect to the mass is 

djCjT^) 0.5 1 _ g 100 e^(T-50) - T-50 

~dhai T00T^^(e™ -1) ~ (e™-i)3 ' ^ ' 

Calculations based on Eq. ([6]) are plotted in Fig. [9], with the vertical axis in arbitrary 

units. 

The ratio of the Einstein oscillator temperature Te to that of the peak in Fig. [9] is 
100/16=6.25 and is independent of Tg. This ratio agrees rather well with the value obtained 
from Figs. 7 and 2 (6.25) for the derivative vs. Mpb and for that vs. Ms from Figs. [8] and 
[2] (6.1). It is interesting to compare these ratios with those obtained for two rather similar 



materials: ZnO (Ref. Il9l ) and GaN (Ref. 



• ij. 



For ZnO one obtains the ratio 7.3 for the 



Mzn derivative and 5.1 for the Mq derivative, whereas in the case of GaN one finds 6.5 for 



the Mca derivative and 3.5 for the Mn derivative. 

We have aheady mentioned in Sec. II that for monatomic crystals the mass derivative of 
the specific heat can be obtained from the corresponding temperature derivative by using 
Eq. (jl])^''. This equation cannot be used for crystals containing several different atoms per 
primitive cell such as PbS, a fact which becomes obvious when one considers that there 
is only one temperature derivative but two (or more) independent mass derivatives. It 



was shown, however, in Ref. 
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(for GaN) and Ref. Il9| (for ZnO) that the sum of the 
corresponding cation and anion mass derivatives equals the result obtained with Eq.(jl]) 
from the temperature derivative. Figure [TO] demonstrates that the same property applies 
to PbS. This figure includes the results obtained by replacing into Eq.(jl]) the experimental 
as well as the ab initio calculated temperature derivatives. The figure also includes the 
two appropriate mass derivatives and their sum (all from the ab initio calculations). Some 
difference appears between the latter and the curve obtained from the measured temperature 
derivative near the maximum at about 20K. We believe that their origin lies in experimental 
and computational errors at these rather low temperatures and also in the neglect of spin- 
orbit coupling in the ab initio calculations. 



VII. CONCLUSIONS 



We have investigated experimentally and theoretically {ab initio) the specific heat of 
PbS and its derivatives with respect to the masses of the two constituents. Reasonable 
agreement is found between experiment and calculations, quantitative differences having 
been attributed to the absence of spin-orbit interaction in the latter. The mass derivatives 
vs. T exhibit maxima at about 12 K for lead and 33 K for sulfur. The large difference in 
the masses of both constituents allows us to separate these maxima and to relate them to 
corresponding maxima in the PDOS. The frequency of the latter is about six times larger 
than that of the former: a simple model, based on a single Einstein oscillator, has been 
proposed to account for this factor. Qualitatively similar ratios have been found in recent 
years for ZnO and GaN, two materials for which the ratio of cation to anion mass is not as 
large as for PbS. Note that the Einstein oscillator curve of Fig|9]can be used, by adjusting 
the oscillator frequency and the vertical scale, to give a rather good representation of the 
experimental data (and the calculations) of Figs. [7| and [HI 
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By inverting the integral equation relating the specific heat and its derivatives to the 
components of the phonon eigenvectors on the two (or more) constituent atoms, it should 
be possible to obtain the corresponding projected densities of states from experimental 
data. Because of the large error bars in some of our measured points, however, we have not 
implemented this procedure, ab initio calculations of the projected PDOS have led us to the 
conclusion, that below 120 cm~^, the eigenvectors are to a large degree lead -like, whereas 
above this frequency they are almost completely sulfur-like. 

In order to remove remaining discrepancies between theory and ab initio calculations 
it should be interesting to modify the ABINIT code so as to include spin-orbit coupling 
in the phonon calculations. Spin-orbit coupling is only included, in the ABINIT version 
available to us for calculating the specific heat, for monatomic materials such as bismuth^i2E. 
A procedure to include spin-orbit splitting in lattice dynamical calculations with ABINIT 
would be highly desirable in order to ascertain the contributions in the phonon dispersion and 
thermodynamical properties in lead chalcogenides, among other multinary semiconductors. 
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FIG. 1: (color online) Phonon dispersion relations of PbS assuming the natural isotope composi- 
tions for Pb and S obtained from the ab initio electronic band structure using the ABINIT code 
(black solid line). Experimental points obtained by inelastic neutron scattering measurements are 
denoted by (blue) diamonds^^. Note that the TO band is considerably higher than the experimen- 
tal data points around the X and K points, a fact which is attributed in the text to the lack of 
spin-orbit interaction in the ABINIT code used^^. 




Frequency (cm"') 

FIG. 2: (color online) Phonon density of states (PDOS) of PbS (total, black solid line) obtained 
from the dispersion relations of Fig. [T]and its components projected on the Pb and S atoms (red 
dots and blue dashes). The horizontal arrows indicate the regions of of predominant Pb and S 
contributions, respectively. For frequencies below (above) 120 cm~^ the sulfur (lead) component 
is extremely small (less then 1%). 
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FIG. 3: (color online) Experimental molar specific heat of PbS (galena, natural isotope composition 
of Pb and S, respectively (red circles). Also shown is the specific heat calculated from the PDOS 
of Fig. [2] (blue solid line) and the data by Parkinson (galena, Ref . J|) . 
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FIG. 4: (color online) Temperature dependence of C^jT'^ as obtained from the ABINIT calculation 
(lower (black) solid curve) and measured by Lykov and Chernik (Ref. 0) (upper dashed (blue) 
curve), by Parkinson (Ref. 8) (black squares), and by us (red circles). The difference between the 
measurements and the ABINIT calculations is attributed to the lack of spin-orbit interaction in 
the latter. 
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FIG. 5: (color online) Effective Debye temperature Od[T) (see text) of PbS. Our data for PbS 
(natural isotope composition for Pb and S, respectively, are displayed by (red) circles. Parkinson's 
data (Ref. |^) by (black) squares, the data reported by Lykov and Chernik (Ref. 9) are represented 
by the (black) solid line. The (green) solid curve represents the results of semiempirical calculations 
by Elcombe (Ref. ,11). 
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FIG. 6: (color online) ABINIT calculations of Cy/T^ for PbS with various isotope compositions, 
(a) "^^Pb'^S with (from top to bottom) n=36, 34, 33, nat, 32 a.m.u. (note that the curves for 
natp^natg ^^^^ natpj-|32g q£ each other). The superscript nat refers to the natural isotope 

abbundance of S which corresponds to an atomic mass of 32.06 a.m.u. (b) ™Pb"^*S with (from 
top to bottom) m=208, nat, 207, 206, 204 a.m.u. The superscript nat refers to the natural isotope 
abbundance of Pb which corresponds to an atomic mass of 207.21 a.m.u. The isotope effects are 
much smaller in (a) than in (b), as corresponds to the low sulfur content of the lower frequency 
eigenvectors. 
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FIG. 7: (color online) ACv/T^ difference for two samples with natural sulfur and two different lead 
isotopic compositions (natural and Mpb = 208 a.m.u.). The (blue) circles represent experimental 
points, the (black) solid curve calculations which are not reliable below ~ 6K. The (red) dotted 
curve represents an extrapolation of the calculated curve for T ^ according to Eq. [5j 




FIG. 8: (color online) AC-u/T^ difference for two samples with natural Pb and two different sulfur 
isotopic compositions (natural and Ms = 34 a.m.u.). The (blue) circles represent experimental 
points, the (black) solid curve calculations which are not reliable below ~ 6K. The (red) dotted 
curve is an extrapolation of the calculated curve for T ^ using Eq. [5j 
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FIG. 9: Einstein oscillator model for the isotope effects of Figs. [7] and El The curve was obtained 
with Eq. ([6]) which assumes, for simplicity, an Einstein temperature of lOOK. The temperature of 
the peak is 16. 5K. The ratio of the Einstein to the peak temperature is, ~6, independent of the 
Einstein temperature. 
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FIG. 10: (color online) The (blue) circles represent the derivative with respect to temperature, 
dln{Cv/T'^)/dlnT of the specific heat of of PbS (galena). The data points were smoothed by 
averaging five neighbor data points on either side (adjacent averages). The (green) solid curve shows 
the r.h.s. of Eq. the (magenta) dashed and (red) dotted curves the derivatives with respect to 
the mass, dln{C^/T^)/dlnM, with M being the mass of either lead or sulphur, respectively. The 
(black) dashed-dotted solid curve is the sum of these derivatives (l.h.s. of Eq. ([1])). Note that 
above 40K the (green) solid and the (black) dashed-dotted curves lie on top of each other. 
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